/** Fichier implémentant une forme de cylindre infini.
 *
 * Ce fichier implémente la classe BaseCylindre.
 *
 * @file BaseCylindre.cpp
 * @author Johan "Solo" Luisier
 *
 * @date 4 Juillet 2011
 */
#include "BaseCylindre.hpp"

#include "Algebre/ConstruitMatrice.hpp"
#include "Algebre/Solveur.hpp"

#include "Utils/StatusCode.hpp"
#include "Utils/debugstreams.hpp"

#include <cmath>
#include <algorithm>
#include <map>

using namespace std;

/** Constructeur par défaut.
 *
 * Ce constructeur construit un cylindre de rayon et hauteur 1, dont
 * "l'origine" est en (0, 0, 0).
 */
BaseCylindre::BaseCylindre()
    : FormeInfinie( "BaseCylindre" ), Axe(), Sommet(), RayonCylindre( 1.0 )
{
    MoteurResolution = new Solveur();
}

/** Constructeur standard.
 *
 * Constructeur permettant la définition de tous les attributs d'une instance
 * de BaseCylindre.
 *
 * @param[in] v axe du cylindre.
 * @param[in] a centre du cercle formant la base de l'isntance.
 * @param[in] r rayon du BaseCylindre.
 */
BaseCylindre::BaseCylindre( const Vecteur& v, const Coord& a, const double& r )
    : FormeInfinie( "BaseCylindre" ), Axe( v ), Sommet( a ), RayonCylindre( abs(r) )
{
    MoteurResolution = new Solveur();

    MatriceSO3 test;

    StatusCode sc( axeZ( Axe, test ) );

    if ( !sc )
    {
        error() << sc.message() << endreq;
        exit( 1 );
    }
    else
    {
        VersGlobal = Transformation( test, Sommet, 1., false );
        VersLocal  = VersGlobal.inverse();
    }
    
}

/** Destructeur, virtuel pour assurer une séquence de destruction correcte.
 *
 * Le destructeur ne fait rien de spécial.
 */
BaseCylindre::~BaseCylindre()
{}

/** Constructeur de copie, permet de copier une instance de BaseCylindre
 * existante.
 *
 * @param[in] cyl BaseCylindre que l'on veut copier.
 */
BaseCylindre::BaseCylindre(const BaseCylindre & cyl)
    : FormeInfinie( "BaseCylindre", cyl.VersGlobal ), Axe( cyl.Axe ),
      Sommet( cyl.Sommet ), RayonCylindre( cyl.RayonCylindre )
{
    MoteurResolution = new Solveur();
}

/** Calcule la valeur du champ de normale au point donné.
 *
 * La normale à un point donné est simplement "l'anti projection orthogonale"
 * du rayon vecteur sur l'axe du cylindre.
 *
 * @warning Ne vérifie pas si le point est bien à la surface.
 *
 * @param[in] pt le point pour lequel on veut la normale.
 *
 * @return la normale au point @c pt.
 */
VecteurU BaseCylindre::normale(Coord const & pt) const
{
    Vecteur tmp(Sommet, pt);

    return VecteurU(tmp.OrthoProject(Axe));
}

/** Permet de vérifier si un point est à @e l'intérieur de l'instance.
 *
 * On calcule simplement la distance entre le point donné et l'axe du cylindre
 * : si cette distance est plus petite ou égale au rayon, alors le point est à
 * l'intérieur du volume.
 *
 * Le calcul est en fait effectué avec le carré de la norme.
 *
 * @param[in] pt point à tester.
 *
 * @retval true si le point est bien dans le volume, @retval false sinon.
 */
bool BaseCylindre::estDans(Coord const & pt) const
{
    Vecteur tmp1(Sommet, pt);
    Vecteur tmp2(tmp1.OrthoProject(Axe));
    return ( tmp2.norme()<RayonCylindre*RayonCylindre );
}

/** Permet de vérifier si un point est à la @e surface de l'instance.
 *
 * On calcule la distance entre le point donné et l'axe, si cette distance est
 * égale au rayon, alors le point est à la surface.
 *
 * Le calcul est en fait effectué avec le carré de la norme.
 *
 * @param[in] pt point à tester.
 *
 * @retval true si le point est bien sur la surface, @retval false sinon.
 */
bool BaseCylindre::estSur(Coord const & pt) const
{
    Vecteur tmp1(Sommet, pt);
    Vecteur tmp2(tmp1.OrthoProject(Axe));
    if ( abs( tmp2.norme() - RayonCylindre*RayonCylindre ) < 1e-7  )
        return true;
    else
        return false;
}

/** Méthode permettant d'appliquer une texture sur le cylindre.
 *
 * On récupère les coordonnées intrinsèques, que l'on traite pour déterminer le
 * pixel à mapper.
 *
 * Si le point en question ne peut être converti en coordonnées intrinsèques,
 * alors le pixel (0,0) est utilisé.
 *
 * @param[in] pt point de la surface à mapper.
 * @param[out] i coordonnée x du pixel de l'image à plaquer au point pt.
 * @param[out] j coordonnée y du pixel de l'image à plaquer au point pt.
 * @param[in] frequence fréquence de répétition de la texture.
 * @param[in] Hauteur hauteur (en pixels) de l'image à mapper.
 * @param[in] Largeur largeur (en pixels) de l'image à mapper.
 */
void BaseCylindre::mapping(const Coord &pt, int &i, int &j,
                           const int &frequence,
                           const int &Hauteur, const int &Largeur) const
{
    Coord carte( intrinseque( pt ) );

    if ( carte == NullePart )
    {
        i = j = 0;
        return;
    }

    double t;
    t = double( frequence * Hauteur ) * carte.x() / ( carte.z() );
    i = static_cast< int >( floor( t ) );

    t = double( frequence * Largeur ) * carte.y() / ( 2. * RayonCylindre );
    j = static_cast< int >( floor( t ) );
}

/** Accessereur pour Axe.
 *
 * Permet de récupérer l'axe du cylindre.
 *
 * @return une référence constante sur Axe.
 */
const VecteurU& BaseCylindre::axe() const
{
    return Axe;
}

/** Accesseur pour Sommet.
 *
 * Permet de récupérer le point d'origine du cylindre.
 *
 * @return une référence constante sur Sommet.
 */
const Coord& BaseCylindre::sommet() const
{
    return Sommet;
}

/** Accesseur pour RayonCylindre.
 *
 * Permet de récupérer la valeur du rayon du cylindre.
 *
 * @return une référence constante sur RayonCylindre.
 */
const double& BaseCylindre::rayonCylindre() const
{
    return RayonCylindre;
}

/** Méthode affichant les paramètres de l'instance.
 *
 * Cette méthode affiche Sommet, Axe, RayonCylindre.
 */
ostream& BaseCylindre::affiche(ostream& os) const
{
    os << "Forme : BaseCylindre" << endl << "Paramètres :" << endl
       << "\tPoint   : " << Sommet << endl
       << "\tAxe     : " << Axe << endl
       << "\tRayon   : " << RayonCylindre << endl;

    return os;
}

/** Calcule l'intersection entre le rayon et le cylindre infini formé par
 * l'instance.
 *
 * On travaille avec les coordonnées locales.
 *
 * Si le rayon est parallèle à l'axe du cylindre, il n'y a pas d'intersection.
 *
 * On résoud l'équation dist( X, a ) = R, avec a la droite définie par l'axe du
 * cylindre et X un élément de la droite passant par le point position avec un
 * vecteur directeur donné par la direction du rayon.
 *
 * A la fin, on a une équation @f$ A t^2 + B t + C @f$.
 *
 * @param[in] ray instance de rayon avec lequel il faut trouver l'intersection.
 *
 * @retval NullePart si l'intersection n'existe pas, @retval inter la
 * coordonnée du point d'intersection.
 */
Coord BaseCylindre::trouveIntersection(const Rayon & ray) const
{
    VecteurU dir( VersLocal( ray.direction() ) );
    Coord orig( VersLocal( ray.origine() ) );

    if ( ( dir ^ Vecteur( 0., 0., 1.) ).norme() < 1e-7 )
    {
        return NullePart;
    }

    double A, B, C;
    double a1( 0. ), a2( 0. ), a3( 1. );
    double d1( dir.x() ), d2( dir.y() ), d3( dir.z() );

    A = a1*a1*( d2*d2 + d3*d3 ) + a2*a2*( d3*d3 + d1*d1 ) +
            a3*a3*( d1*d1 + d2*d2 ) - 2 * ( a1*a2*d1*d2 + a1*a3*d1*d3 + a3*a2*d3*d2 );

    B = 2 * ( a1*a1 * ( d2 * orig.y() + d3 * orig.z() ) +
              a2*a2 * ( d3 * orig.z() + d1 * orig.x() ) +
              a3*a3 * ( d1 * orig.x() + d2 * orig.y() ) ) -
            2 * ( a1 * a2 * ( d1 * orig.y() + d2 * orig.x() ) +
                  a1 * a3 * ( d1 * orig.z() + d3 * orig.x() ) +
                  a2 * a3 * ( d2 * orig.z() + d3 * orig.y() ) );

    C = a1 * a1 * ( orig.y() * orig.y() + orig.z() * orig.z() ) +
            a2 * a2 * ( orig.z() * orig.z() + orig.x() * orig.x() ) +
            a3 * a3 * ( orig.x() * orig.x() + orig.y() * orig.y() ) -
            2 * ( a1 * a2 * orig.x() * orig.y() +
                  a1 * a3 * orig.x() * orig.z() +
                  a3 * a2 * orig.z() * orig.y() ) - RayonCylindre * RayonCylindre;

    MoteurResolution -> setParams( C, B, A );

    vector< double > solutions;

    MoteurResolution -> resoud( solutions );

    vector< double >::const_iterator solIt, solFin( solutions.end() );

    Vecteur dist;

    for (solIt = solutions.begin(); solIt != solFin; solIt++)
    {
        if ( *solIt > 1.e-7 )
        {
            Coord pt( orig + (*solIt) * dir );
            dist = Vecteur( pt );
            pt = VersGlobal( pt );
            return pt;
        }
    }

    return NullePart;
}

/** Méthode calculant les coordonnées intrinsèques d'un point de la surface.
 *
 * Si le point n'est pas sur la surface, on renvoie NullePart.
 *
 * On va retourner un angle et une longeur : la longueur est 0 si le point est
 * confondu avec Sommet, elle est égale au rayon pour un point situé sur le
 * corps du cylindre à une hauteur nulle. Si le point est confondu avec le
 * symétrique de Sommet la longueur est 2 * RayonCylindre + mHauteur.
 *
 * @param[in] pt point de la surface.
 *
 * @retval NullePart si le point n'est pas à la surface, @retval local dont les
 * coordonnées sont @code ( theta * r, l, r ) @endcode
 */
Coord BaseCylindre::intrinseque(const Coord& pt) const
{
    if ( ! estSur( pt ) )
        return NullePart;

    double theta, l;

    Vecteur v( VersLocal( pt ) );

    l = v * VecteurU( 0., 0., 1. );

    VecteurU unit( v );

    theta = acos( unit.x() );

    if ( v.y() < -1.e-7 )
        theta = 2 * Forme::ValeurPi - theta;

    return Coord( theta * sqrt( pow( v.x(), 2 ) + pow( v.y(), 2 ) ),
                  l, sqrt( pow( v.x(), 2 ) + pow( v.y(), 2 ) ) );
}
